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We investigate numerically in spherical geometry the interaction of stratification with 
precession. Both stable stratification and unstable stratification are studied. In the pa- 
rameter regime we are concerned with, stable stratification suppresses the precessional 
instability, whereas unstable stratification and precession can either stablise or destablise 
each other at the different precession rates. 
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1. Introduction and motivation 

The magnetic field of many celestial bodies is maintained by the dynamo effect which 
requires that the body contains a fluid conductor executing a non-trivial motion. The en- 
ergy source for this motion is frequently assumed to be thermal convection. Accor dingly, 
convection driven dynamos in spherical spheres have been extensivel y studied (iBussel 
l2000t iRoberts k Glatzmaierl [2005 IWicht k Tilgnerll2010t |Jonedl201ll ). iBullardl (lll) llili 
first pointed out that precession is a possible energy source for the geodynamo. Precession 
driven flows, their stabi lity and dynamo properties have been investigated independently 
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of thermal convection ( Poincar eill910l; lMalkuslll9 68; Buss ej|l968t IVanvo k Likins 
l999otlNoir et ai.ll200lUTilgner k Bussej|200lULorenzani k Tilgnerll200l[ 
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More recently, some attention has also been paid to other types of me- 



chanical forcing, such as tides, libration and collisions (jWeiss et aZ.ll2008h iLe Bars et al 



201lUDwver et q/.ll201lh 



For the Earth and the other planets and moons of the solar system, the mechanical 
forcing is well known, while it is uncertain if their cores are convecting. In general, it will 
always be necessary to consider buoyancy and a mechanical forcing such as precession 
in conjunction. The thermal stratification can be either unstable or stable depending on 
which epoch in the thermal history of a celestial body is being considered. Typically, 
the surface cooling of a recently formed body leads to a superadiabatic gradient. As 
the heat generated during accretion and formation of the body is progressively lost, the 
internal temperature gradient eventually drops below the adiabatic gradient, yielding 
a stable thermal stratification. Small bodies such as planetesimals or the Moon have 
already reached this point in the past. 

In this paper, we investigate the hydrodynamics of stably and unstably stratified fluid 
in a precessing spherical shell. This is intended as a model of the Earth, whose precession 
is maintained by the gravitational torque of the Moon and the Sun, but also as a model 
of planetesimals in the early solar system which undergo force free precessional motion 
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after collisions if the direction of their angular momentum does not coincide with one of 
their principal axes of inertia. 

In section 2 we formulate the problem and introduce the numerical methods. In section 
3 we investigate the interaction of stable stratification and precession. In section 4 we 
investigate the interaction of unstable stratification and precession. We conclude with 
some discussion in section 5. 



2. Equations and numerical methods 

We consider a fluid in a spherical shell of inner radius and outer radius r a . Suppose 
that the spherical shell spins at the angular velocity fl s and precesses at the angular 
velocity £l p . A background temperature TJ, is imposed on the fluid, which we assume to 
be a linear function of radius, 

T b = ^-^(r-r )+T , (2.1) 
a 

where Tj and T a are the temperature at and r Q , respectively, and d — r Q — n is the 
thickness of the gap. Such a background temperature is maintained by a heat source 
or sink which varies inversely proportional to radius (because V 2 T& = 2(T Q — Ti)/(rd)). 
An alternative model in which the stratification is imposed by fixed temperatures at 
the boundaries in the absence of a heat source leads to a non-uniform stratification and 
has not been considered for simplicity. In the frame moving with the boundary, the 
dimensional Navier-Stokes equation of a Boussinesq fluid reads 

— + U ■ VlX + rt re f X (fl re f X 7") + 2fl re f X U + (fl p X fi s ) X T = 

- —Vp + vV 2 u+ —g, Vw = 0, (2.2) 

Po Po 

where Cl re f = J~i s + £l p is the angular velocity of the moving frame with respect to the 
inertial frame and p the fluid density at r . In the Boussinesq approximation, density 
variations are only retained in the buoyancy term and the density deviation is propor- 
tional to the temperature deviation, 

P^ = -a(T-T o ) = -a(0 + T h -T o ), (2.3) 

Po 

where = T — Tf, is the temperature deviation from the conductive profile and a 
the thermal expansion. In a sphere of constant density, the gravitational acceleration 
is proportional to radius, 

9 = -go — , (2.4) 

To 

where g a is the gravity at r a . 

Substituting (|2.3p and (|2.4[) into (|2.2[) and normalising length with d, time with fij 1 , 
u with £l s d and with (Tj — T D ), we obtain the dimensionless Navier-Stokes equation 

du ~ * * — - 

— h u • Vu = -V$ + EkV 2 u + 2u x (z + Pofl p ) + Po(z x fl p ) x r + i?a0r, (2.5) 

where all the curl- free terms, e.g. the centrifugal force and the term associated with Tf,, 
are already absorbed into the total pressure $. In (|2.5|) there are three dimensionless 
parameters. The Ekman number 
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measures the ratio of the viscous time scale to the spin time scale, the Poincare number 

P„=a ,2.7) 

measures the precession rate. Despite its unusual form, we choose to call the parameter 
in the buoyancy term the Rayleigh number 

K _5*g^a, (2. 8) 

which is proportional to the imposed stratification. A negative Ra corresponds to a stable 
stratification and a positive Ra to an unstable stratification. For a stable stratification 
Ra is just the square of the dimensionless Brunt- Vaisala, frequency. The unit vector along 
the precession axis in Cartesian coordinates (x, y, z) is 

ti p = sin /3 cos i a: — sin/3sinty + cos j3z , (2-9) 

where z is the spin axis and /3 the angle between the spin axis z and the precession axis 
ftp. The boundary condition of velocity is no-slip at the outer boundary, whereas it is 
stress-free at the inner boundary to approximate a full sphere. 

Accordingly, the dimensionless equation of temperature deviation is 

£^ + M . V e- Ur = — v 2 e. (2.10) 

at Pr 

In (j2~10|) the Prandtl number 

Pr = - (2.11) 

K 

measures the ratio of viscosity to the thermal diffusivity. The boundary condition is 
9 = at both inner and outer boundaries. 

Equations (|2.5j) and (|2.10|) are solved numerically. In the calculations we fix the aspect 
ratio ri/r a to 0.1 such that the spherical shell is almost a spherical cavity and hence the 
inner core is almost negligible, and the angle /? to 60°. We investigate only retrograde 
precession (Po < 0) which is geophysically relevant. We fix the Prandtl number to 
Pr = 1. In the calculations of stable stratification we fix Po = —0.3 and vary Ek and Ra 
to investigate how the stable stratification influences the precessional instability. In the 
calculations of unstable stratification, we select two Ekman numbers Ek = 5x 10~ 3 and 
1 x 10~ 3 j_which are sufficiently high to maintain the flow precessionally stable, and vary 
Po and Ra to investigate the interaction of precessional and convective instabilities. 

The numerical calculations are carried out in the spherica l coordinates (r, 0) with 



the parallel pseudo-spectral code provided by iTilgnerl (| 1999 6h . The toroidal-poloidal de- 
composition method is employed such that the divergence free condition of fluid flow 
V • u = is automatically satisfied. The spherical harmonics P[ n (cos8)e lm ^ > are used on 
the spherical surface (#,</>) and the Chebyshev polynomials Tfc(r) are used in the radial 
direction. Resolutions as high as 256 (radial r), 128 (colatitude 9) and 64 (longitude (/>) are 
used. A semi-implicit scheme is employed for time stepping, using an Adams-Bashforth 
scheme for the nonlinear and a Crank-Nicolson scheme for the diffusive terms. 



3. Stably stratified precessional flow 

The precessional flow in the interior of spherical container is mainly a solid body 
rotation with angular velocity Clf which is in general different from both f2 s and fi p 
(|Busselir968l ). At the boundary the fluid rotation matches the spin rate fi s such that 
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there exists an Ekman boundary layer, where the poloidal flow is generated by Ekman 
pumping. In addition to the bo undary layers th ere also exist internal shear layers that are 
spawned at the critical latitude ( Tilgnerll 1999a ). A well- understood instability m echanism 
in these flows is a triad resonance between the basic flow and tw o inertial modes (IKerswell 
19951 iTilgner fc BussdlioOll lLorenzani fc Tilgner]|200l[ 



We now consider the interaction of a stable stratification and precession. Radial strat- 
ification tends to suppress fluid motion in the radial direction, which suggests that stable 
stratification suppresses precessional instability. However, stable stratification leads to 
gravito-inertial waves whose frequencies may be closer to a perfect triad resonance than 
the pure inertial waves of th e unstratified c ase, so that destabilisation through stable 
stratification is possible, too ( Kerswell 19931) . Only stabilisation is observed in our nu- 
merical calculations. In thejiumerical calculations we fix Po = —0.3 and calculatejyarious 
combinations of Ek and Ra, namely Ek ranges from 5 x 10~ 3 to 1 x 10~ 4 and Ra from 
(purely precessional flow) to —1. Because a stable precessi onal flow is anti-symmetric 
about the centre, i.e. u{— r) = —u(r), we use symmetries ( Lorenzani fc Tilgner l200"ll) 
to detect the precessional instability, i.e. the flow is separated into two parts, u s (r) = 
[u(r) — u(— r)} /2 and u a (r) — [u(r) + u{— r)] /2, and the kinetic energy of u a repre- 
sents the instability. 

Figure [T] shows the total kinetic energy E (figure 1(a)) and the ratio of instability 
energy E a to the total energy (figure 1(b) ) against the Ekman number Ek at different 
Rayleigh numbers Ra (for time dependent flows, E and E a are calculated by averaging 
over time). Figure 1(a) shows that increasing Ek reduces the energy E. This behavior is 
well known from previous studies and is due to the fact that the viscous term is more 
important at higher Ek which damps any motion excited by the Poincare force. It is 
also seen from figure 1(a) that a larger \Ra\ reduces the energy. This is consistent with 
figure 2(a) which shows the modulus of the fluid rotation vector |u>| = 

to indicate that a lower \Ra\ corresponds to a stronger fluid rotation. In ord er to explain 
the data in figures 1(a) and [21 one would have to extend Busse's calculation ( Busse 19681) 
to a stably stratified medium, which is not attempted here. 

Figure [2] shows the radial dependence of the fluid rotation vector. The shear in the basic 
flow is due to viscous corrections to the solution of the inviscid equation, which is simply 
a solid-body rotation. Viscosity introduces Ekman pumps and shear layers crossing the 
entire fluid volume, all of which possess a radial velocity component. Since stratification 
suppresses radial motion, it also reduces the shear. In addition to the modulus of the 
fluid rotation vector, we also calculated the angle 7 between the fluid rotation vector 
and its average, 7 = arccos(u; • ZJ / (\u>\\uj\) , and 7 is less than around 10° in the interior 



and reaches around 30° near the boundary. Moreover, it is interesting that figure 2(b) 



indicates a counter- rotation of the fluid in the precession frame (lj z + 1 < 0), which is 
located in the vicinity of the inner boundary. 

The orientation of the fluid rotation axis is best characterised in the precession frame in 
which the geographic and precession axes are stationary. Figure |3] shows the angle formed 
by the rotation axis of the fluid (averaged over the fluid volume) with the geographic axis 
(figure 3(a) ) and with the meridian of the precession axis (figure 3(b) ) in the precession 
frame. Figure 3(a) parallels the total kinetic energy shown in figure 1(a) since the kinetic 
energy increases with the increasing angle between the fluid and boundary rotation axes. 
The key point in figure 3(a) is that the stable stratification does not noticeably modify the 
rotation axis until Ra = — 1, and therefore we may conclude that a stable stratification 
takes effect on the precessional flow only for \Ra\ > 1. 
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(a) (b) 

Figure 1. Stably stratified precessional flow at Poincare number Po = —0.3. The total kinetic 
energy E (a) and the ratio of instability energy E a to total energy (b) as a function of Ekman 
number Ek at different Rayleigh numbers Ra. 




(a) (b) 

Figure 2. Stably stratified precessional flow. Radial dependence of fluid rotation vector at 
Ekman number Ek = 2 x 10 -3 and Poincare number Po — —0.3. Modulus \oj\ (a) and z 
component ui z (b). 



The energy in the unstable modes is shown in figure 1 (b) E a ^ indicates instability. 

The onset of instability is shifted towards smaller Ek if \Ra\ is increased. This can be 
seen in more detail in figure 4(a) which shows the critical Ekman number, below which 
the flow is unstable, as a function of Ra. The critical Ekman number monotonously de- 
creases with increasing \Ra\. Even though this is intuitive (because stable stratification 
suppresses radial motion and hence retards the onset of instability), the opposite could 
have happened as well if the changes in frequency of the inertial modes due to strat- 
ification had brought a combination of them closer to a triad resonance. There is an 
indication of this effect at finite amplitude, since the hierarchy among the different Ra 
is not preserved in figure 1(b) at Ek = 1CP 3 . 

The onset of the instability is to a large degree just a matter of the energy in the 
basic flow. The growth rate of a triad resonance depends on the shear in the basic flow 
which increases with the energy of the basic flow. This pointjs demonstrated in figure 
|4(b)| which shows that the onsets of instability for different Ra nearly coincide if E a is 
plotted as a function of E, with the exception of Ra = — 1. This again shows that stable 
stratification is significant only for \Ra\ — 1 or larger. 




Figure 3. Stably stratified precessional flow at Poincare number Po — —0.3. Position of fluid 
rotation axis as a function of Ekman number Ek at different Rayleigh numbers Ra. Colatitude 
(a) and longitude (b) of fluid rotation axis in the precession frame. The dashed line denotes the 
position of precession axis. 




Figure 4. Stably stratified precessional flow at Poincare number Po = —0.3. (a) The squares 
show points determined numerically to lie on the stability limit. The points are connected by 
a line to guide the eye. (b) The instability energy Ea as a function of the total energy E for 

different Rayleigh numbers Ra. 



4. Unstably stratified precessional flow 

For the study of the fluid core of most planets it is likely more relevant to investigate 
the interaction of unstable stratification with precession. In iLe Bars fc Le Dized ( 2006) 
the geometry of an infinitely long cylinder with an elliptical cross section was analytically 
investigated and it was shown that the elliptical instability and the convective instability 
may either stabilise or destabilise each other in different parameter regimes. In this section 
we numerically calculate unstably stratified precessional flow in spherical geometry to 
study the interaction of precession with unstable stratification. 

We select two Ekman numbers, Ek = 5 x 10~ 3 and 1 x 10~ 3 , such that the purely 
precessional flow is stable, namely it is stable at least until \Po\ — 1 for Ek = 5 x 10~ 3 
and \Po\ — 0.3 for Ek = 1 x 10~ 3 . Then we test various combinations of Po and Ra to 
seek the neutral stablity curves of the critical Rayleigh number Ra c against \Po\ at both 
Ekman numbers, as well as to calculate^ some supercritical flows at Ek = 5 x 10~ 3 . Ra c 
is sought by increasing or decreasing Ra in steps of 0.01 at a given Po . We distinguish 
two different types of flows by the dominant azimuthal wavenumbers in the spectra 
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of the velocity field. These spectra are computed in a frame of reference in which the 
axes of precession and boundary rotation are stationary and a z'-axis is pointing along 
the rotation axis of the fluid. The azimuthal angle <p' is measured around the z'-axis 
and the azimuthal wavenumbers are denoted by m ■'. The onset of convection occurs at 
m' = 3, so that we call the flow convective if the kinetic energy is concentrated at 
mf = 3. If on the contrary the spectrum has its largest contributions at m! < 3, we call 
the flow precessional. Flows undergoing precessional instability have their kinetic energy 
concentrated in these wavenumbers. The motion with wavenumbers ml = 1 and 2 can be 



identi fied as two inertial modes by the same technique as used by lLorenzani fc Tilgner 
( 20011 ). These two modes, together with the spin-over mode, fulfill the requirements of 
a triad resonance. Their parameters (latitudinal wavenumber V and frequency zu') are 
close to two analytically det ermined eigenmod es with m! = 1, V = 6, w' = —0.537 and 
m' = 2,1' = 6, w' = -1.093 (| Greenspan! 1 196i ). 



Figure 5(a) shows an approximative stability diagram of unstably stratified preces- 
sional flows at Ek = 5 x 10 -3 . The solid line denotes the neutral stability curve above/below 
which flow is unstable/stable, and the circle symbol denotes the convective flows and the 
square symbol the precessional flows. In the absence of precession (Po = 0), Ra c for 
the onset of convective instability is 0.45 and the dominant mode is m! = 3. When \Po\ 
increases to 0.1, Ra c decreases and the instability is still convective. The instability be- 
comes precessional for \Po\ = 0.2 or larger. Ra c decreases until it reaches a minimum of 
Ra c = 0.19 at \Po\ = 0.45. The flow becomes more stable with stronger precession for 
\Po\ > 0.45. At \Po\ = 0.6 the flow is so stable that Ra c reaches 1. At such a high Ra 
the flow pattern is convective again. 

Precession has apparently a dual role. On one hand, both precession and convection 
can lead to instabilities on their own so that a superposition of both can be expected to 
be even less stable. This is observed in figure 5(a) around —Po ~ 0.45, where the critical 



Rayleigh number is reduced by more than a factor of 2 by precession, but precession 
alone in the isothermal fluid (corresponding to Ra = 0) is stable. On the other hand, 
precession introduces shear on top of the global rotation, and the onset of convective 
instability has now to be computed for a sheared basic flow, which can lead to a higher 
critical Rayleigh number. An example of this phenomenon is seen in figure 5(a) for —Po 
around 0.6. 

We then investigate the lower Ekman number of Ek = 1 x 10~ 3 . We only seek the 
neutral stability curve but do not calculate any supercritical flows, and we calculate until 
\Po\ = 0.3 beyond which the purely precessional flow (Ra = 0) is already unstable (figure 



1(b)). Figure 5(b) shows this neutral stability curve. As in figure 5(a) the circle symbol 
denotes the convective instability and the square symbol the precessional instability. For 
purely convective flows (Po = 0), the critical Rayleigh number Ra c = 0.07 at Ek = 
1 x 10~ 3 is much lower than Ra r = 0.45 at Ek = 5 x 10~ 3 . This is consistent with the 



asym ptotic scaling law of previou s calculation s (Roberts 1968; Busse 



20071 ) and numerical calc ulations (IZhang et aLll2007t iTilgner fc Busse 



197C : IZhang et al 



19971 ). Notice that 



the definition of Ra in ( Tilgner fc Bussd 19971) should be translated to our definition 



(equation 12. 8|) such that the asymptotic scaling law is Ra c = 0(Ek 2 / 3 ) (this scaling law 
does not precisely hold in our numerical calculations because we choose too high Ekman 
numbers). For Ek = 1 x 10~ 3 , Ra c has a maximum value of Ra c = 0.21 at \Po\ = 0.2, 
which shows that precession stabilises the flow at small \Po\ whereas it destabilises at 
large \Po\. We know that the neutral stability curve at Ek — 5 x 10~ 3 (figure 5(a)) 



eventually decreases when \Po\ is large enough for precession alone to be unstable. In 




Figure 5. Diagram of convective stability and precessional stability for unstably stratified flow 
at Ek — 5 x f0~ 3 (a) and Ek = 1 X fO -3 (b). Points on the solid line denotes points on the 
neutral stability curve. The circle symbols denote the convective flows and the square symbols 
the precessional flows. 



summary, at both the high and low Ekman numbers the neutral stability curve is not 
monotonic but it has a minimum at the high Ek and a maximum at the low Ek. 

To end this section we discuss the Nusselt number Nu which measures the ratio of 
the total heat transfer to the thermal conduction in the fluid at rest. The stratifying 
linear temperature profile assumed here is maintained by a heat source so that the Nus- 
selt number depends on radius. The Nusselt number is conveniently computed at the 
boundaries as 

^=^« (4.1) 

The brackets denote an average over the spherical surface of radius r Q or and time. 
Because of the internal heat source, the Nusselt numbers at the two boundaries are 
related through 

Nu\ r = ro - rfNu\ r=rt = 1 - f? 2 , (4.2) 

where ry = r,i/r — 0.1 is the aspect ratio used in our calculations. It is verified by 
our numerical calculations that equation (14.21) precisely holds. Since Nu at the inner 
boundary can be deduced from the one at the outer boundary, we only show Nu at the 
outer boundary. Figure |5] shows Nu at the outer boundary against Po at Ek — 5 x 10 -3 
for Ra = and Ra = 0.5. Because a purely precessional flow (Ra — 0) has some poloidal 
component it also transfers heat, in which case the temperature deviation is a passive 
scalar. Figure [6] contains several points with Nu around 1.05 which are compared in table 
[T] At equal poloidal kinetic energy, one flow may be more efficient at tranporting heat 
than the other because of a better correlation between radial velocity and temperature. 
According to this criterion, convection is more efficient than precession at advecting heat. 
Adding precession to a convective flow reduces this correlation by a factor of about 4, as 
shown in table [1] 



5. Discussion 

We investigated the interaction of precession with thermal stratification. The heat 
transport in precessional flow is less than in convective flow at equal rms of the radial 
velocity component due to a smaller correlation between radial velocity and tempera- 
ture. Unstable stratification together with precession can be either more stable or more 
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Figure 6. Nusselt number Nu at the outer boundary as a function of Po at Ekman number 
Ek = 5 x 10~ 3 and Rayleigh numbers Ra = and 0.5. 
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Table 1. Correlations (J u r OdV), {J u^dV) and (J Q 2 dV) of purely precessional and convective 
flows. The brackets denote the time average. The first four flows at Ra = are precessionally 
stable and the next three flows at Ra = 0.5 are convectively unstable as shown in figure [5(a)] 



unstable than the stratification or precession acting alone. The same could have been ex- 
pected from the combination of precession with stable stratification, but in this case, we 
only found examples in which the stratification stabilises the flow. The presence of stable 
stratification becomes relevant for precession if the ratio of the Brunt- Vaisala, frequency 
to the rotation frequency is near 1. 

We present a rough estimate to determine what this criterion implies for various celes- 
tial bodies. Consider the extreme case of a body which has cooled down so much that it 
is nearly isothermal. Within the Boussinesq approximation used here, this corresponds 
to a stable stratification with a temperature gradient equal to the adiabatic gradient, 
given by 

=-^. (") 

. 9z J ad C P 

where g is the gravitational acceleration, a the thermal expansion coefficient and c p the 
specific heat capacity at constant pressure. The gravitational accleration at the surface 
of a sphere with radius R is given by 

g = ^TtGpR, (5.2) 



where G is the gravitational constant. The Brunt- Vaisala frequency for the cooled body 
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under consideration is given by 

N = \r^Wh = ga {i = r GpaR {l- (5 - 3) 

If we accept the following values to be representative of planetary cores 

p = 11390 kg/m 3 , a = 1.5 x lO^K" 1 , T = 4000K, c p — 860 J/ (kg • K), (5.4) 
we obtain 

N = 1.03 x Hr^-RnrV 1 . (5.5) 

If we take as an example planetesimals with a typical radius of R = 10 km or R = 100 km, 
we find N — 10 _6 s _1 or N = 10~ 5 s ~ 1 , respectively. T he rotation of planetesimals 
has a typical period of about 10 hours dWeiss et al. 2008), i.e. the rotation rate is fi 



1.7 x 10 -4 s _1 , which implies that the buoyancy is negligible in this context. The situation 
is different for the Earth with a core of radius 3.4 x 10 6 m and a rotation period of 24 hours. 
A stabilising stratification with a gradient 4% less than the adiabatic gradient would be 
enough for the stratification to significantly influence the flow driven by precession. 
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